N.est=function(data,n,t,s,estmethod){
  data=as.name(data); n=n; t=t; s=s
  estmethod=as.name(estmethod)
  if(estmethod=='direct'){
    est=ceiling(n*t/s)
    vhat=(t^2*n*(n-s))/s^3
    bound=round(2*sqrt(vhat),2)
    lower=ceiling(est-bound); upper=ceiling(est+bound)
  } else if(estmethod=='Chapman'){
    est=ceiling(((t+1)*(n+1)/(s+1))-1)
    vhat=((t+1)*(n+1)*(t-s)*(n-s))/((s+1)^2*(s+2))
    bound=round(2*sqrt(vhat),2)
    lower=ceiling(est-bound); upper=ceiling(est+bound)
  } else if(estmethod=='inverse'){
    est=ceiling(n*t/s)
    vhat=(t^2*n*(n-s))/(s^2*(s+1))
    bound=round(2*sqrt(vhat),2)
    lower=ceiling(est-bound); upper=ceiling(est+bound)
  }
  cat("","\n","Population estimation results from method:",estmethod,'\n',"Data:",data,"\n","n =",n,'t =',t,'s =',s,"\n",
      "Estimate of N =",est,"\n","Vhat =",vhat,"\n","Bound =",bound,"\n","Lower Bound =",lower,"Upper Bound =",upper,"\n","")
  results=list(estmethod=estmethod,data=data,n=n,t=t,s=s,est=est,vhat=vhat,bound=bound,lower=lower,upper=upper)
}
# to use the function with its call:
# N.est(data,n,t,s,estmethod)

# data: name of dataset, in quotes
# n: sample size, single value
# t: number of tagged items, single value
# s: number tagged in 2nd sample, single value 
# estmethod: c('direct','Chapman','inverse')